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Abstract 

The anisotropy of temperature is studied here in a strong two-dimensional shockwave, sim- 
ulated with conventional molecular dynamics. Several forms of the kinetic temperature are 
considered, corresponding to different choices for the local instantaneous stream velocity. A 
local particle-based definition omitting any "self" contribution to the stream velocity gives the 
best results. The configurational temperature is not useful for this shockwave problem. Config- 
urational temperature is subject to a shear instability and can give local negative temperatures 
in the vicinity of the shock front. The decay of sinusoidal shockfront perturbations shows that 
strong two-dimensional planar Shockwaves are stable to such perturbations. 
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I. INTRODUCTION 



Shockwaves are useful tools for the understanding of material behavior far from 
equilibrium^! 2 -^. The high-pressure physicist Percy Bridgman played a key role in the 
adaptation of experimental shockwave physics to the thermodynamic characterization 
of materials at high pressure^. Because Shockwaves join two purely equilibrium states, 
shown to the left and right of the central shockwave in Figure 1, the experimental and 
computational difficulties associated with imposing nonequilibrium boundary conditions 
are absent. 

Begin by assuming that the flow is both stationary and one- dimensional. Such a flow 
gives conservation of the mass, momentum, and energy fluxes throughout the shockwave. 
Thus the fluxes of mass, momentum, and energy, 

pu, P xx + pu 2 , pu[e + (P xx /p) + (u 2 /2)] + Q x , 

are constant throughout the system, even in the far-from-equilibrium states within the 
shockwave. Here u is the flow velocity, the velocity in the x direction, the direction of 
propagation. We use conventional notation here, p for density, P for the pressure tensor, e 
for the internal energy per unit mass, and Q x for the component of heat flux in the shock 
direction, x. It is important to recognize that both the pressure tensor and the heat flux 
vector are measured in a local coordinate frame moving with the local fluid velocity u(x). 

In a different "comoving frame" , this time moving with the shock velocity and centered 
on the shockwave, cold material enters from the left, with speed u s (the "shock speed") 
and hot material exits at the right, with speed u s — u p , where u p is the "piston speed". 
Computer simulations of Shockwaves, using molecular dynamics, have a history of more 
than 50 years, dating back to the development of fast computers^&^&i&ii. Increasingly 
sophisticated high-pressure shockwave experiments have been carried out since the Second 
World Wari 2 -. 

In laboratory experiments it is convenient to measure the two speeds, u s and u p . These 
two values, together with the initial "cold" values of the density, pressure, and energy, 
make it possible to solve the three conservation equations for the "hot values" of p, P xx , 
and e. A linear relation between P xx (x) and V(x) = 1/p, results when the mass flux 
M = pu is substituted into the equation for momentum conservation: 

P xx (x) + M 2 /p(x) = P hot + (M 2 /p hot ) = P cold + (M 2 /p cold ) . 

2 



Figure 1 











Us 


Us - Up 



Figure 1: A stationary one-dimensional shockwave. Cold material enters at the left with a speed 
u s , passes through the shockfront which separates the cold material from the hot, and exits at 
the right with speed u s — u p . The cold-to-hot conversion process is irreversible and corresponds 
to an overall entropy increase. 

This nonequilibrium pressure-tensor relation is called the "Rayleigh Line". The energy 
conservation relation, 

Ae = -(l/2)[P hot + P cold ]Ay , 

based on equating the work of compression to the gain in internal energy, is called the 
"Shock Hugoniot Relation". See Figure 2 for both. A recent comprehensive review of 
shockwave physics can be found in Ref. 3. 

Laboratory experiments based on this approach have detailed the equations of state for 
many materials. Pressures in excess of 6TPa (sixty megabars) have been characterized^. 
If a constitutive relation is assumed for the nonequilibrium anisotropic parts of the pres- 
sure tensor and heat flux, the conservation relations give ordinary differential equations 
for the shockwave profiles. With Newtonian viscosity and Fourier heat conduction the 
resulting "Navier- Stokes" profiles have shockwidths on the order of the mean free path^ 8 -. 

Early theoretical analyses of Shockwaves emphasized solutions of the Boltzmann equa- 
tion. Mott-Smith's approximate solution of that equation- 1 —, based on the weighted aver- 
age of two equilibrium Maxwellian distributions, one hot and one cold, revealed a tem- 
perature maximum (and a corresponding entropy maximum) at the shock center, for 
shocks with a Mach number exceeding two. The Mach number is the ratio of the shock 
speed to the sound speed. Twenty years later, molecular dynamics simulations showed 
shockwidths of just a few atomic diameters.- 1 ^* 8 - These narrow Shockwaves agreed nicely 
with the predictions of the Navier-Stokes equations for relatively weak shocks. At higher 
pressures there is a tendency for the Navier-Stokes profiles to underestimate the shock- 
width. Some of the computer simulations have shown the temperature maximum at the 
shock front predicted by Mott-Smith 8,10 . These temperature maxima are constitutive 
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Figure 2 
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Figure 2: Calculated Rayleigh Line and Shock Hugoniot relation for the simple repulsive van 
der Waals' equation described in the text. The Rayleigh line includes nonequilibrium states 
within the shock while the Shock Hugoniot line is the locus of all equilibrium states accessible 
by shocking the initial state. 

embarrassments — they imply a negative heat conductivity for a part of the shockwave 
profile. 

Our interest here is twofold. We want to check on the stability of planar shock- 
waves (the foregoing analysis assumes this stability) and we also want to characterize 
the anisotropy of temperature in the shock. This latter topic is particularly interest- 
ing now in view of the several definitions of temperature applied to molecular dynam- 
ics simulations^^^^£^. A configurational-temperature definition as well as several 
kinetic-temperature definitions can all be applied to the shockwave problem. 

The plan of the present work is as follows. Section II describes the material model 
and computational setup of the simulations. Section III describes the computation of the 
shock profiles and the analysis confirming their stability. Section IV compares the various 
kinetic and configurational temperature definitions for a strong stable shockwave. Section 
V details our conclusions. 
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II. SHOCKWAVE SIMULATION FOR A SIMPLE MODEL SYSTEM 



We simulate the geometry of Figure 1 by introducing equally-spaced columns of cold 
particles from a square lattice. The initial lattice moves to the right at speed u s . The 
interior of the system is purely Newtonian, without any boundary, constraint, or driving 
forces. Those particles coming within the range of the forces (er = 1) of the righthand 
boundary have their velocities set equal to u s — u p (the mean exit velocity) and are 
discarded once they reach the boundary. This boundary condition, because it corresponds 
to a diffusive heat sink at the exit, only affects the flow in the vicinity of that boundary. 
We initially chose the speeds u s = 2 and u p — 1 to correspond to twofold compression. 
The approximate thermal and mechanical equations of state, 

e=(p/2)+T; P = pe . 

together with the initial values, 

Pcoid = 1 ; T co i d = ; P co i d = (1/2) ; e cold = (1/2) , 

give corresponding "hot" values. These two sets of thermodynamic data mutually satisfy 
the Rayleigh Line and Shock Hugoniot Relation for twofold compression from the cold 

state: 

p hot = 2 ; T hot = (1/4) ; P hot = (5/2) ; e hot = (5/4) . 

The mass, momentum, and energy fluxes are 2, (9/2), and 6, respectively. For reasons 
explained below it was necessary to modify these conditions slightly, using instead u s = 
2u p = 1.75. 

We considered a wide variety of system lengths and widths and found no significant 
difference in the nature of the results. For convenience we show here results for a system 
of length L x = 200 and width L y = 40. Because the density increases by a factor of two 
in the center of the system the number of particles used in this case is about 12000. The 
number varies during the simulation as new particles enter and old ones are discarded. 
The total length of the run is one shock traversal time, 200/m s , though the shock is 
itself localized near the center of the system and in fact moves very little in our chosen 
coordinate frame. In retrospect, the simulations could just as well have been carried out 
with a much smaller L x . Because we wished to study stability we felt it necessary to use 
a relatively wide system. 
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For a two-dimensional classical model with a weak van der Waals' repulsion the equa- 
tion of state described above follows from the simple canonical partition function: 

Z l/N oc VTe~ p/2T ; PV/NkT = [d\nZ/d\nV} T ; E/NkT = [dhxZ/d\nT] v ■ 

Our original intent was to use the smooth repulsive pair potential 

<p = (10/trx 2 )[1 - (r/a)] 3 for r < a — > 

($) ~ (Np/2) / (j)(r)2Trrdr = Np/2 . 
Jv 

for a sufficiently large a (3 or so) that the simple equation of state was accurate. But for 
large a the shockwidth is also so large that detailed studies are impractical. In the end we 
chose to set the range of the forces equal to unity, so that the initial pressure is actually 
zero rather than 1/2. Nevertheless, the choice of u s = 2 is still roughly compatible with 
twofold compression. Equilibrium molecular dynamics simulations for a = 1 give the 
following solution to the conservation relations for twofold compression, from po = 1 to 
p = 2 with u s = 2u p = 1.75: 

pu : 1.0 x 1.75 = 2.0 x 0.875 = 1.75 ; 

P + pu 2 : 0.0 + 1.0 x 1.75 2 = 1.531 + 2.0 x 0.875 2 = 3.062 ; 

pu[e+(P/p) + (u 2 /2)} = 1.75[0.0 + 0.0 + 1.531] = 1.75(0.383 + 0.766 + 0.383] = 1.75 x 1.531 . 

We introduced a sinusoidal perturbation in the initial conditions by using twice as 
many columns of particles (per unit length) to the right of a line near the center of the 
system: 

^shock = 6sm(2iry/L y ) . 

The time development of a system starting with this sinewave displacement perturbation 
is shown in Figure 3. The panel corresponding to the time t = Dt shows that the decay 
is underdamped, while the following panels show that the planar shockwave is stable. 

Figure 4 illustrates the effect of the gradual underdamped flattening of the sinewave 
perturbation on the density profile. The propagation of the wave is followed through five 
shock traversal times of the observation window used in Figure 3. 

The density profiles are computed with Lucy's one-dimensional weighting function 2 ^! 

w lD {r < h) = (5/4/i)[l - 6(r/h) 2 + 8(r//i) 3 - 3(r//i) 4 ] ; r < h = 3 
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Figure 3: Particle positions in the initial condition correspond to t = 0. Those to the left of 
the sinewave boundary move to the right at speed u s = 1.75 while those to the right travel at 
speed u s — u p = u s /2. Particle positions at five later equally-spaced times are shown too. The 
time interval Dt = 40 /u s is 2000 timesteps. The 40 x 40 windows shown here would contain 
1600 cold particles or 3200 hot particles at the cold and hot densities of 1.0 and 2.0. 
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Figure 4: Density profiles at the same times as those illustrated in Figure 3. These one- 
dimensional density profiles were computed with Lucy's one-dimensional smooth-particle weight 
function using h = 3 and the particle coordinates shown in Figure 3. Increasing time corresponds 
to increasing line thickness as the Shockwave moves slowly to the left. 
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Figure 5 
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Figure 5: Density profiles exactly as in Figure 4, but with a reduced range, h = 2. The 
shockwave moves slowly to the left. Note the wiggly structure, a consequence of choosing h too 
small.. 



for all combinations of Particle i and gridpoint k separated by no more than h = 3. 
For comparison density profiles using a shorter range, h = 2, are shown in Figure 5. 
These latter profiles exhibit a wiggly structure indicating a deterioration of the averaging 
process. 

Similar conclusions, using similar techniques, have been drawn by Hardy and his 
coworkers*^, who were evidently unaware of Lucy and Monaghan's work. In their inter- 
esting analysis of a two-dimensional shockwave, Root, Hardy, and Swanson use a weight 
function like Lucy's, but with square (rather than circular) symmetry and with only a 
single vanishing derivative at its maximum ranged They discuss the ambiguities of de- 
termining temperature away from equilibrium and rightly conclude that the range of the 
spatial weight function needs careful consideration. Hardy- recently told us that his use 
of a spatial weighting function was motivated by a conversation with Philippe Choquard. 

The profiles we show in Figure 4 are typical, and show that the shockwidth rapidly 
attains a value of about 3 particle diameters, and has no further tendency to change 




The density is computed by evaluating the expression 



; w lk = wl (\xi-x k \) , 
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as time goes on. A detailed study shows that the sinewave amplitude exhibits under- 
damped oscillations on its way to planarity. Evidently, for this two-dimensional problem, 
the one-dimensional shockwave structure is stable. In the next Section we consider the 
temperature profiles for the stationary shockwave. 

III. CONFIGURATIONAL AND KINETIC TEMPERATURES IN THE 
SHOCKWAVE 

The kinetic and configurational contributions to temperature and pressure have been 
discussed and explored in a variety of nonequilibrium contexts. Shockwaves, with sta- 
tionary boundary conditions far from the shockfront, allow the anisotropy of the kinetic 
temperature to be explored, analyzed, and characterized with purely Newtonian molecu- 
lar dynamics. Kinetic-theory temperature is based on the notion of an equilibrium ideal 
gas thermometer— i22j2i. The temperature(s) measured by such a thermometer are given 
by the second moments of the velocity distribution: 

{kT x K x ,kT y K y } = m{(vl),(v 2 y )} 

where the velocities are measured in the comoving frame, the frame moving at the mean 
velocity of the fluid. 

A dilute Maxwell-Boltzmann gas of small hard parallel cubes undergoing impulsive 
collisions with a large test particle provides an explicit model tensor thermometer for the 
test-particle kinetic temperature^. The main difficulty associated with kinetic tempera- 
ture lies in estimating the mean stream velocity, with respect to which thermal fluctuations 
define kinetic temperature. In what follows we compare several such definitions, in an 
effort to identify the best approach. 

Configurational temperature is more complicated and lacks a definite microscopic me- 
chanical model of a thermometer able to measure it. Configurational temperature is based 
on linking two canonical-ensemble equilibrium averages, as was written down by Landau 
and Lifshitz more than 50 years ago^: 

7i is the Hamiltonian governing particle motion. Unlike kinetic temperature this con- 
figurational definition is independent of stream velocity. Its apparent dependence on 
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Figure 6: Particles with negative configurational temperatures are shown here, using the data 
underlying Figs. 4 and 5. The material to the right of the shock is "hot", with cold material 
entering at the left. The smaller dots correspond to negative values of T® x and the larger ones 
to negative Tyy. 

rotation^ is negligibly small for the (irrotational) shockwave problems considered here. 
We apply both the kinetic and the configurational approaches to temperature measure- 
ment here, considering individual degrees of freedom within a nominally one-dimensional 
shockwave in two-dimensional plane geometry. 

Although the notion of configurational temperature can be defended at equilibrium, 
the shockwave problem indicates a serious deficiency in the concept. An isolated row or 
column of particles, all with the same y or x coordinate and interacting with repulsive 
forces is clearly unstable to transverse perturbations. The symptom of this instability is 
a negative "force constant" 



or yy 



< . 



The uncertain sign of VV7i explains the presence of wild fluctuations, and even negative 
configurational temperatures, in the vicinity of the shockfront. See Fig. 6. This outlandish 
behavior means that the configurational temperature is not a useful concept for such 
problems. 

Defining the kinetic temperature requires first of all an average velocity, about which 
the thermal fluctuations can be computed. The velocity average can be a one-dimensional 



10 



sum over particles {i} sufficiently close to the gridpoint k: 

i i 

Alternatively a local velocity can be denned at the location of each particle i by using a 
two-dimensional Lucy's weight function, 

w 2D (r <h) = (5/tt/i 2 )[1 - 6{r/h) 2 + 8{r/h) 3 - 3(r//i) 4 ] ; r < h = 3 . 

h 

27rrw(r)dr = 1 , 
and summing over nearby particles {j}: 

j 3 

These latter two-body sums can either both include or both omit the "self" term with 
i = j. Our numerical results support the intuition that it is best to omit both the self 
terms^. We compute all three of these x and y temperature sets for the stationary 
shockwave profile and plot the results in Figure 7 and 8. 

The data show that in every case T xx exceeds T yy in the leading edge of the shock 
front. There is a brief time lag between the leading rise of the longitudinal temperature 
T xx and the consequent rise of the transverse temperature T yy . This is to be expected 
from the nature of the shock process, which converts momentum in the x direction into 
heatj2&2Z The Rayleigh line itself (see again Figure 2) shows the mechanical analog of 
this anisotropy, with P xx greatly exceeding P yy . Note that the lower pressure in Figure 2, 
the Hugoniot pressure, is a set of equilibrium values, (P xx + P yy )/2. It is noteworthy that 
the height of the temperature maximum is sensitive to the definition of the local velocity. 
Evidently a temperature based on the local velocity, near the particle in question, and 
with a coarser averaging range h = 3, gives a smaller gradient and should accordingly 
provide a much simpler modeling challenge. 

IV. SUMMARY 

This work shows that planar Shockwaves are stable for a smooth repulsive potential 
in a dense fluid. We also find that a mechanical instability makes the configurational 
temperature quite useless for such problems. Our investigation of temperature definitions 

11 




Figure 7 
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Figure 7: Typical instantaneous temperature profiles, at t = 5Dt for the strong shockwave 
described in the text. Local particle-based definitions of stream velocity give the two somewhat 
lower temperature pairs, (T^ X ,T^}. The grid-based definitions of T^ x and T^ y are indicated 
with the heaviest lines. They both use a one-dimensional weight function. This definition gives 
a strong temperature maximum for T*~,. In all three cases the longitudinal temperature T^ x 
exceeds the transverse temperature T^ y near the shock front. The temperatures obtained with 
two-dimensional weights and including the "self" terms in the average velocity are indicated by 
the dashed lines and are significantly lower than the rest throughout the hot fluid exiting the 
shockwave. The "correct" hot temperature, far from the shock, is kT = 0.13, based on separate 
equilibrium molecular dynamics simulations. The profiles shown here were all computed from 
the instantaneous state of the system after a simulation time of 5Dt = 200 /u s . 

shows that relatively smooth instantaneous temperature profiles can be based on the 
weight functions used in smooth particle applied mechanics. It is noteworthy that the 
constitutive relations describing inhomogeneous nonequilibrium systems must necessarily 
include an averaging recipe for the constitutive properties. The present work supports 
the idea^ that the range of smooth averages should be at least a few particle diameters. 

The better behavior of a particle-based temperature when the "self" terms are left out 
is not magic. Consider the equilibrium case of a motionless fluid at kinetic temperature 
T. The temperature of a particle in such a system should be measured in a motionless 
frame. If the "self" velocity is included in determining the frame velocity an unnecessary 
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Figure 8 
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Figure 8: Temperature profiles with the same data as in Figure 7, but with a reduced averaging 
range, h = 2. Again, the temperatures including "self" contributions to the local velocity are 
shown as dashed lines. Notice the sensitivity of the maximum in T^ x to the range h. Far to the 
right, the "correct" hot temperature, at equilibrium, away from the shock, is kT = 0.13, based 
on separate equilibrium molecular dynamics simulations. 

error will occur. Thus it is plausible that the "self" terms should be left out. The data 
shown in Figs. 7 and 8 support this view. 

The instantaneous particle-based velocity average provides a smoother gentler profile 
which should be simpler to model. By using an elliptical weight function, much wider in 
the y direction than the x, one could consider the grid-based temperature as a limiting 
case. The somewhat smoother behavior of the particle-based temperature recommends 
against taking this limit. The elliptical weight function would leave the divergent config- 
urational temperature unchanged. 
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